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Abstract. The Vlasov equation is well known to provide a good description of the dynamics 
of mean-field systems in the N —> oo limit. This equation has an infinity of stationary states 
and the case of homogeneous states, for which the single-particle distribution function is 
independent of the spatial variable, is well characterized analytically. On the other hand, 
the inhomogeneous case often requires some approximations for an analytical treatment: 
the dynamics is then best treated in action-angle variables, and the potential generating 
inhomogeneity is generally very complex in these new variables. We here treat analytically 
the linear stability of toy-models where the inhomogeneity is created by an external field. 
Transforming the Vlasov equation into action-angle variables, we derive a dispersion relation 
that we accomplish to solve for both the growth rate of the instability and the stability threshold 
for two specific models: the Hamiltonian Mean-Field model with additional asymmetry and 
the mean-field 4 model. The results are compared with numerical simulations of the A^-body 
dynamics. When the inhomogeneous state is stationary stable, we expect to observe in the 
Af-body dynamics Quasi-Stationary-States (QSS), whose lifetime diverges algebraically with 
N. 
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1. Introduction 

Long-range forces can be found in a wide variety of physical systems, including self- 
gravitating systems, Coulomb systems, wave-plasma interactions and two-dimensional 
hydrodynamics. The interest in studying long-range forces has been revived in the last 
decade, not only because of the broad domain of physical systems involving such forces, 
but also because of the presence of unusual phenomena, both at equilibrium and out of 
equilibrium. Let us mention negative specific heat, temperature jumps, broken ergodicity 
and quasi-stationary states. Reviews and books have been recently published in this field 

E El 12 011319. 

A particular, but interesting, case is the one of mean-field interactions, for which each 
particle is directly coupled to all the others with equal strength, whatever their distance. 
Although this is an idealization, it serves as a useful approximation and appears, in addition, 
to give at least the good trend. Moreover, there are physical situations in which particles are all 
in interaction via a field, whose dynamics is in turn determined by the motion of the particles 
themselves: this is for example the case of wave-particle interactions in plasmas 0, Free 
Electron Lasers 0, Collective Atomic Recoil Lasers 10 and Traveling Wave Tubes iTTOll . 
This self-consistent effect can also be obtained in systems composed only of particles by 
introducing a coupling to an order parameter, as it is done for the Hamiltonian Mean Field 
(HMF) model lfTTl[T2l[T3ll . which has been widely studied in recent years as a paradigm for 
systems with long-range interactions [Q]| . 

The kinetics of models with N particles and only mean-field interactions is exactly 
described, in the infinite N limit, by the Vlasov equation |[T4l [L51 . This equation exhibits 
an infinity of stationary solutions and its dynamical evolution starting from a generic initial 
state can be extremely complex. Focusing on stationary states, their stability has been studied 
using different methods, but mainly by restricting the analysis to homogeneous stationary 
states, that are characterized by a single-particle distribution function which is independent 
of the spatial variable. These states are of major interest in kinetic theory, because they often 
constitute the "supposed" physical equilibrium state. For instance a globally neutral plasma 
has an equilibrium which is also locally neutral, giving a homogeneous charge distribution. If 
perturbed, this state is expected to be stable, showing a relaxation back to the homogeneous 
state ruled by Landau damping [fT6l [P71 [T8Tl . This phenomenology is also observed in the 
HMF model |fT9l , for which the homogeneous state is stable above a given energy threshold, 
which depends on the initial momentum distribution. 

However, below this energy, the homogeneous state is unstable and one observes a 
dynamical evolution towards inhomogeneous states, whose stability properties are much more 
difficult to determine. Inhomogeneous states appear for example in gravitational dynamics 
ll20l . because of the attractive nature of the Newton force. Their stability has been studied in 
the context of the Vlasov equation, yet the necessity to resort to action-angle variables j[2T| 
makes the problem analytically tricky. Apart from numerical approaches (see e.g. E2l ). 
one can project the dynamics onto a Fourier basis, yet at a cost of performing infinite 
sums 11251 l26l : then, only a truncation can yield tractable results. Such technique was also 
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used in the context of plasmas [1231 12311 . where the waves often generate inhomogeneous 
states; expanding the dynamics along modes, such as Hermite polynomials BTTl . requires 
anyhow a truncation in the sums. Analytical results were also obtained on BGK modes, 
whose stability properties were connected, in the small inhomogeneity limit, to those of 
homogeneous states 11281 129~1 . Later on, the unstable nature of periodic BGK modes under 
specific perturbations was rigorously shown OOllBTI . but the problem remains open for other 
types of systems and perturbations. More recently, some general criteria were proposed to 
derive the stability of inhomogeneous states 0211331 . 

Some toy models were also studied whose states are naturally inhomogeneous: this is 
typically the case of systems when an external potential is present in addition to the self- 
consistent one OS [351 EH EH. A first interesting model 03 E21 ES EH is the mean-field 
ip 4 model: an Ising-like spin variable is represented by a scalar field in one dimension, acted 
upon externally by a double- well potential which selects two states; the mean-field term of the 
Hamiltonian is a quadratic coupling of the scalar field at two different lattice sites. A second 
interesting model is a generalized version of the Hamiltonian Mean-Field (HMF) model to 
which an anisotropic external potential is added 0711 . 

In this paper we focus the above mentioned toy models, and show that one can treat 
exactly the stability of inhomogeneous states. The Vlasov equation will be rewritten in action- 
angle variables [1391140114111 and we will focus on those inhomogeneous stationary states whose 
single-particle distribution function does not depend on the angle variable, i.e. those that are 
homogeneous in angle. We will derive a general stability criterion which, besides giving the 
value of the threshold energy (action) at which these stationary states destabilize, will allow 
us to obtain the growth-rate of the instability. 

In Section[2jwe will introduce and discuss the Vlasov equation in action-angle variables 
and we will derive the stability condition for inhomogeneous states and for generic mean-field 
and external potentials. In Sections E] and |4j we shall apply the general method introduced in 
Section[2jto the specific cases of the anisotropic HMF model and of the mean-field </? 4 model, 
deriving explicit analytical expressions for the stability threshold and for the growth rate of the 
instability. These theoretical predictions will be then compared with numerical simulations 
performed with iV-body Hamiltonians. Finally, in Section [51 we will draw some conclusions 
and we will discuss some perspectives of this work. 

2. The Vlasov equation in action-angle variables and the stability relations 

Let us consider N particles in one-dimension whose positions and momenta are (qj,pj), 
j = l,...,N. They interact through the two-body (symmetric) potential v(qj,qk) and, in 
addition, each particle is trapped into the external potential W(qj). Hamilton's equations for 
such a system are 



4j = Ph 

pj= -W(qj)-d qj V[{q k )](qj), 



(1) 

(2) 
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where V [{qt}] {qj) = U/AO v(qj, qu) stands for the mean-field potential acting on particle 
j. The 1 /N term is a rescaling factor Il42ll which allows one to perform the mean-field limit 
discussed in Refs. |[l4~l[T5l . The prime will denote, from now on, the derivative with respect 
to the position variable q. Eqs. (OQ) and © can be derived from the following Hamiltonian 



V, i A 



(3) 



ff = 2y + %) + /W](?i) 
j ^ 

where the (qj,Pj) are couples of canonically conjugated variables. Let us introduce the so- 
called empirical measure 

1 N 

f(9, P^"> = ^Yj 6{q ~ 9jiO)S(p - Pj (t)) . (4) 

7-1 

It can be shown |[T5l that, in the N — » oo limit, the single-particle distribution function 
f(q, p, t) obeys the following Vlasov equation 

drf + pd q f - (W'(q) + V'[f](q)) d p f = , (5) 

where 

V[f](q, t) = ff dq'dp' f(q', p' , t)v(q, q'), (6) 

is the averaged mean-field potential. One can also show that the yV-body dynamics is well 
described by the Vlasov equation over times that are at least of order lnA^ lfl5l . This makes 
the Vlasov framework a natural one to study such systems when a large number of particles 
is involved. 

The Vlasov equation can also be written in Hamiltonian form using the following 
functional 

H[f] = jf dqdp f(q, p, t) UL + W(q) + 1 V[f](q)\ . (7) 

After having introduced the appropriate Poisson brackets for the functionals A[f] and B[f] 

' d 5A d SB d 8A d 6B 
K dp8fdq6f dq 8f dp Sf 



r A m rC J J ^ l d 5A d 5B d 5A d 8B\ 
{A,B} = JJ dqdp f(q,p,t)[———— -————] 



the dynamics of A[f] is given by 

d t A = {H,A}. (9) 

If one rewrites the single particle distribution function in the functional form f(q, p, t) = 
ff dq'dp' f(q', p', t)6(q - q')S(p - p'), one obtains the evolution equation 

d,f(q,p,t) + - 

dp dq dq dp 

= d t f(q, p, t) + {h[f](q, p), f(q, p, t)} = (10) 

where h[f](q, p) = p 2 /2 + W(q) + V[f](q) and the brackets are now the standard Poisson 
brackets. This equation is nothing but the Vlasov equation ©. 

It is straightforward to check that the Boltzmann-Gibbs equilibrium distribution 
/bg(<7>P) = Z -1 exp(-8h(q, p)), with 6 an arbitrary constant and Z a normalization constant, 



Stability of inhomogeneous states 



5 



is a stationary solution of this equation (i.e. d t f B c = 0). In fact, all distributions that depend 
on (q, p) only through h are stationary. The existence of an infinity of stationary distributions 
is actually responsible for the peculiar out-of-equilibrium regimes in which Af-body long- 
range systems get trapped over very long times 0]|. More specifically, starting from a generic 
unstable distribution, a long-range system typically relaxes towards a "quasi- stationary" state, 
which can be significantly different from Boltzmann-Gibbs equilibrium. Quasi-Stationary 
States (QSS) can be interpreted as stable stationary states of the Vlasov equation in the N —* oo 
limit. The relaxation to statistical equilibrium occurs on much longer time scales, that were 
observed to diverge either algebraically |[T9l or logarithmically 0711 with N, depending on 
whether the "quasi-stationary" state corresponds to a stable or an unstable stationary state of 
the Vlasov equation. Relaxation to equilibrium is not due to collisions but due to finite-iV 
effects (also called "granularity"), which can be modeled by convenient kinetic equations, 
like Landau or Lenard-Balescu equations 03 Q/71 dU. Stable stationary states of the Vlasov 
equation are therefore of paramount importance in order to understand the dynamics of long- 
range systems. It is therefore crucial to determine the general conditions for stationarity and 
stability, for both homogeneous and inhomogeneous states. 

Let us consider the stationary state fo(q, p). If one focuses on the Lagrangian trajectory 
of a single particle, one immediately realizes that it is a constant energy trajectory of the 
energy functional 

P 2 

h[f ](q,p) = y + W(q) + V[f ](q), (11) 

which is a straightforward consequence of Eqs. (TJJ) and ©. Hence, it is convenient to cast the 
dynamics into the appropriate variables associated with this trajectory, namely the "action- 
angle" variables 

J(h) = ^ j) p(h, q')dq' = i- j) ^2(h-W(q')-V[f ](q'))dq' (12) 

dq' 

yj2(h-W(q')-V[fo](q')) 
where the frequency to is given by 



(/> =<o \ - Uq (13) 

Jo 



OJ 



i r d 4 



dh 

dl' 



(14) 



2n J ^2(h-W(q')-V[f ](q')), 
It is important to note that the conjugate variables (/, (p) are not action-angle stricto 
sensu: Since Vlasov dynamics is infinite dimensional and only a specific set of conserved 
quantities can be typically identified (e.g. the Hamiltonian, total momentum, the Casimirs 
ff dqdp C(f(q,p)), with C an analytic function), its integrability is not generic |[43l . The 
term action-angle variables comes from the fact that the dynamics of a Lagrangian test- 
particle is integrable if the single-particle distribution function is stationary. Indeed, for a 
stationary distribution f , the potential V[f ] is constant in time. Therefore, the dynamics 
of the test-particle is that of a one-degree-of-freedom system with the associated conserved 
quantity h[f ], hence integrable. A dependence of the potential on time caused by a non- 
stationary distribution f(q, p, t) would introduce an extra 1 /2 degree of freedom, thus breaking 
integrability a priori. 
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In this single particle framework and for stationary distributions, the energy h depends 
only on the action J, so that a particle evolves on a trajectory of constant "action" J at 
the constant action-dependent angular speed (p = djh(J) = co(J). The change of variables 
(q, p) —* ((f), J) being canonical, the corresponding Poisson brackets, which apply to functions 
of the phase-space, are equivalent 

{a, b} q , p = d p ad q b - d q ad p b = {a, b)^j = djad^b - d^adjb. (15) 

Using this equivalence and the condition d$h = 0, the Vlasov equation (0) for f can be 
recast in the following form 

dMWtfo = a>(J)dtf = . (16) 

Hence, the stationarity condition, d t f = 0, leads to f = fo(J)- This means in particular that 
the stationary distributions are those that are homogeneous in angle, with any distribution in 
action J. Such a result highlights the relevance of action-angle variables for the analysis of 
Vlasov stationary dynamics, but also for the study of QSS. 

We shall now consider a perturbation 6f around f , that is f(cp, J) = f (J) + 5f(cf), J). The 
linearity of the potential V with respect to the distribution, as emphasized by its definition in 
Eq. ©, implies that V[f] = V[f ] + V[Sf]. Using property (UJ) for the Vlasov equation (151101) 
and neglecting second-order terms in 6f leads to the linearized Vlasov equation 

d t 6f + o*J)df6f - (d p f )V'[6f](<f>, J) = 0, (17) 

where the factor d p f should be expressed in terms of (<p, J) and the derivative of V is with 
respect to q and then it is also expressed in terms of (<p,J). The study of this equation in 
full generality would imply the solution of an initial value problem using a Laplace-Fourier 
transform and then a transformation back to action-angle variables using a Bromwich contour 
[fT71 [TSl . We will be here less ambitious and we will focus on the study of an eigenmode 
8f{(j),J;t) = e At f((p,J) with the eigenvalue A determining the stability properties. Inserting 
this ansatz solution in Eq. (fTTl) one gets 

(A + a>(/)<W - (d p fo)V [f] (cp, J) = 0. (18) 

Assuming a non-zero co (the frequency a> typically only vanishes on the separatrices of the 
single particle phase-space), the above equation turns into 

^ ( e W)/) _ __ (5p/o)r [y] (0> 7) = o. (19) 

After integration over the angle tp, and assuming that the integration constant vanishes, one 
gets 

-A4>/o(.J) rt> 

f 77- tye#M(d p f W,J)V'\f](<l>',J) = 0. (20) 

oj(J) Jo L J 

This equation can be fully cast into action-angle variables using the following relation 
df 0f , dJdMJ) dhdJ p 

-£-(.<!> P) = 74 TT~ = ^~^fo( J ) = -fo(J) , (21) 

op op oJ op oh 10 
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f = -trAfV^ (25) 



which, inserted into Eq. results in the following dispersion relation 

f-fo(J)— 277T #W,^)^ Ml7) V' / (f,J) = 0. (22) 
^> 2 (^) Jo L J 

It is convenient to express the integral in this latter equation in terms of the position variable 

q' . Indeed, using Eq. (PT3T ), the differential d(p' can be calculated as a function of q' at constant 

action J, which means along a single-particle trajectory. One gets 

d(f> = = = —dq , (23) 

yl2(h-W(q')-V[f ](q')) P 

which allows one to put Eq. (1221) into the following form 

e -A4>/a)(J) rq 

f - foi-D—nr eA *' MJ) Y \A teW = °> < 24 ) 

U{J) Jo L J 

in which the integral is performed at constant action J. The interest of this alternative formula 
is that it may be easier to solve in some cases. In particular, if one focuses on the stability 
threshold, given by taking A = 0, the integral over q' can be solved straightforwardly and 
Eq. (1241) can be rewritten as 

co(J) 

where q is, in general, a function of both action and angle. 

Since all functions in angle are 2^-periodic, it is common to project the dispersion 
relation in a Fourier base ll3~9ll4TT| . However, since in Eq. (|22l) both the p term and the potential 
V[f] have generically a non trivial dependence on the angles, one ends up with expressions 
where all Fourier modes are coupled. The modes are decoupled only when momentum does 
not depend on angle, which is the case of homogeneous states, for which momentum coincide 
with action (modulo a sign). 

In what follows we will discuss a method which allows us to compute the stability 
threshold and the growth rate A without resorting to a Fourier expansion. The method is, 
however, not generic and its application depends on the specific form of the interaction 
potential. We will therefore discuss separately two examples. 

3. The HMF model with additional asymmetry 

Introduced in Ref. lT3"71l . the HMF model with additional cosine on-site potential is a 
generalization of the paradigmatic HMF model [|TTL [T2l [T3l . Besides the mean-field term 
v(qj, q k ) = - cos (qj - q k ), an external potential W of amplitude k is present 

W(qj) = k cos 2 qj . (26) 

The Hamiltonian © reads 

H[f] = JJ dqdpf(q,p) y ■ + /ccos 2 ^-^ jj dq'dp' f(q',p') cos (q - q' 

At variance with the HMF model, the spatially-homogeneous state is no longer a stationary 
state of the Vlasov equation, due to the presence of the on-site potential. 



') 



.(27) 
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Using formula (12"2~1) . one easily gets the dispersion relation for this model 



where 



e A<f,/aj(J) pp _ _ 

f-foW—iFK W pW> J) ^ MJ) \ M *W sin ^'' J)) ~ M yW cos (^ ■ / )) = o, (28) 
u 2 (J) Jo L 

M[/]=M,[/] + /M,[/] = ^J* dqdpf(q,p)cosq + iJf dqdpf(q, p) sin ^(29) 

stands for the magnetization. For the sake of simplicity, and p' will respectively refer 
to q((p, J), q((f>', J) and p((f>', J) in the remaining of this section. Equation (1281) can be solved 
by multiplying each term by either cos q or sin q, and then integrating over phase-space. One 
gets the following equations 

M x [f] (l - IxAM) + M y [f]I A xx [f ] = 0, (30) 

-M x [f]I* Y [fo] + M y [f] (l + I* x [f ]) = 0, (31) 

where 

Ixj[fo] = J dJ^j)d<pe-W^X(q)f q dq'e A *'^Y(q'), (32) 

and the label X (resp. Y) stands for the cos (resp. sin) function. The integral j> is performed 
over a single-particle trajectory. 

Inhomogeneous stationary states of the Vlasov equation correspond to solutions of the 
linear system of equations (130H31I) with non vanishing (M x , M y ). They can be found only 
when the determinant vanishes. This condition allows to rewrite the dispersion relation in the 
form 

(l - Ix, Y m) (l + Iixlfo]) + IyjlfoVixifo] = 0. (33) 

The numerical resolution of this equation can be performed by using the explicit expressions 
of the action-angle coordinates P31 , for a particle of energy h and position q 



Jinih) = 



n 



Jout(h) = & 7 <Pout(q,h) = (35) 

n \h) 2 Jt {Kin) 

where the label in I out stands for inside/outside of the separatrix of the potential k cos 2 q, while 
$ , J£T and & are elliptic integrals of the first kind. 

In order to compute the growth rate Re(A) from Eq. (1331) it is necessary to choose a 
specific unperturbed stationary distribution fo(J). We here consider "waterbag" distributions 
in action-angle space that are homogeneous in angle: these are two-level distributions, which 
are nonzero and homogeneous between two lines of constant action J = Ji and J = J 2 

2n(J 2 - Ji) 

where the first factor guarantees the normalization of the density f , while is the Heaviside 
step function. Morever, we here focus on waterbags delimited by a given energy U, i.e. we 
consider all trajectories with energies h < U (so that J 2 = J 2 (U) and = 0), such as those 
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Figure 1. Waterbags in action-angle (panel a) and in (q, p) space (panel b). The waterbags 
have increasing boundary energies U — 0.2, 0.4 and 0.55 and they are represented by filled 
contours of lighter and lighter grey as the energy is increased. The dashed line corresponds to 
the separatrix, which has energy U s = 0.3 and action J s = 0.5. 



0.8 




0.2 0.3 0.4 0.5 0.6 

U 



Figure 2. Growth rate Re(A) (full line) of the instability of the inhomogeneous waterbag states 
obtained by solving Eq. ( 1331 for waterbags with boundary energy U. The crosses are the results 
of numerical simulations of the A^-body Hamiltonian. The agreement between theory (which 
decribes the N — > oo limit) and numerics (which is performed at N = 3 x 10 5 ) is reasonably 
good apart from the region near the separatrix energy U s = 0.3 and the one near the critical 
energy U c = 0.498, which is theoretically determined by solving Eq. ( l42l . 



represented in Fig. Ufa) and (b). It is interesting to remark that, since the change of variables 
(q, p) «-» (0, J) is canonical, f (q, p) is also a two-step distribution with the boundary given 
by the curve h(q, p) = U. It should be pointed out that, altough the action fixes the energy 
univocally, a trajectory of given energy is always splitted in two: those with positive and 
negative momentum p for U > U s = 0.3, the separatrix energy, and the ones with < q < n 
and n < q < In for U < U s . This has the consequence that, when performing integrations 
over the action-angle space, the two trajectories give separate contributions. Related to this 
remark is for example the evaluation of the normalization of fa. the total area of the waterbag 
is indeed 2 x 2n(J(U) - 7(0)) = AnJ(U). 

For the waterbag initial conditions, the integral in Eq. (1321) reads 

1 



2n{J 2 -J x ) 



— £d(P e~^ luj{Jl) X{q) f 7 dq' e^' MJl) Y(q') 

^(■A) J Jo 
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— Sdcpe-^'^Xiq) f dq'e^' MJ2) Y(q') 
(J2) J Jo 



•(37) 

The numerical solution of Eq. (1331) . using Eq. (1371) . are then compared with the result 
of simulations performed with the Af-body Hamiltonian using a sixth-order integration 
scheme ll46l with time step 0.1. Figure [2] shows the growth rate Re(/l) obtained theoretically 
(full line) as a function of the boundary energy U . The growth rate is determined numerically 
by fitting an exponential to the short-time increase of the magnetization. One notices the 
existence of a threshold energy U c = 0.498 (determined more precisely in the following), 
which separates a region where the waterbag is stable (U > U c ) from one where the waterbag 
is unstable (U < U c , Re(/t) > 0). When the waterbag is stable, the iV-body dynamics shows 
a QSS regime with zero magnetization but with an inhomogeneous distribution of particles 
in the q spatial coordinate. Let us remark that the theoretical results shows a divergence of 
Re(/l) at the separatrix energy U = U s = 0.3 where the frequency a)(J s ) = 0: this divergence 
is not reproduced by the Af-body dynamics. Moreover, in the iV-body dynamics, the threshold 
energy is found to be around U ~ 0.44, well below the theoretical value. Indeed, in the energy 
region 0.44 < U < U c the growth of the magnetization is spoiled by finite-iV effects, and its 
exponential character is not clear any more. However, the energy U c is really the one where 
we numerically observe a de stabilization of the zero magnetization state. 

The critical energy U c beyond which the waterbags become stable can be explicitly 
derived using Eq. (1331) and by imposing A = 0. Let us first note that, in this equation, the last 
term vanishes, since both I xx [fo] and 7 yy [/o] yield an integral of sing cos g over a trajectory. 
Consequently, the product (1 - /£ y [/o])(l + I^xifo]) should be zero. Then, considering that 
\p\ = V2£/ c -\/l - (k/U c ) cos 2 q, integrating over q', and using Eq. (1231) and then Eq. (PT21) . we 
finally get 

I°xAfo] = j^j £dcf>cos 2 q (38) 



1 f , cos 2 q 

K dq (39) 



2ttJ, 



cos 2 q 



\IU C -K cos 2 q 

(40) 



2 j> \dq\ yjU c - kcos 2 q 



1 fm- 



sin 2 q 



n V U c -K cos 2 q 

JU/b] = - o7 / = ■ ( 41 ) 

1 <p \dq\ -y U c — k cos 2 q 

Let us explain the meaning of the uncommon notation \dq\. When integrating over segments 
of the single-particle trajectory where p is negative, q decreases. Thus, both dq and p are 
negative, so that their ratio or product is positive. The use of the differential \dq\ allows us to 
unify notation for both the cases in which p and dq are positive or negative. The coefficient 2 
in front of the first integral originates from the double boundary of the waterbag, be it inside 
or outside the separatrix. It can be shown that both expressions (1401) and (|4T|) are strictly 
decreasing functions of U c . Moreover, integral (1401) tends to one in the U c — » k limit, so that 
1 - I XY [fo] is always positive. The threshold of stability is thus given by solving the implicit 
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equation 



sin 2 g 



= 2 \dq\ yJU c -k cos 2 q. 



(42) 



a/L^ - KCOS 2 q 

The numerical resolution of the above equation for k = 0.3 yields the value U c ~ 0.498, in 
excellent agreement with the energy value at which Re(/t) vanishes (see Fig.©. 

We note that the above derivation of the threshold energy U c corroborates with the result 
derived in Ref. 071l . where the same result was obtained by developing the single-particle 
distribution as a sum of derivatives of Dirac distributions. The truncation of the expansion to 
the very first term allowed the authors of Ref. Il37l to obtain the same implicit equation (142)) . 
The approach presented here is more general, since it provides a dispersion relation for any 
stationary distribution, and allows us to derive the stability condition without any additional 
hypothesis. 

We devote the final part of this Section to the derivation of the growth rate of the 
instability and of the threshold energy for the HMF model, in the limit where the on-site 
potential is turned off (k = 0). Although this result was already obtained [0111441 . its derivation 
in this new context allows us to point out the connection between action-angle variables 
(4>, J) and the canonical ones (q,p). In fact, when only the mean-field potential couples the 
particles, the non-magnetized inhomogenerous stationary states become homogeneous in q 
and, in correspondence, the action-angle variables reduce, modulo a sign, to the canonical 
coordinates 



J = (j) P(h, q) dq 
dh 

° J= dJ = H 



\P\, 



C q dq' 
J P 



sign(p) q. 



(43) 
(44) 
(45) 



The presence of absolute values is due to the fact that the action-angle variables take into 
account the direction of the motion along the trajectories, which are now ballistic. Then, 
inserting the following relations 



f 



sing'dg' = 



,Aqlp 



1 + A 2 /p 2 



A . 
— smq 

P 



cos q 



M ' lp 



cos q'dq' 



e' 



MIp 



A 

sin q H — cos q 
P 



1 + A 2 /p 2 

into Eq. (1321) . one can explicitely write the dispersion relation (1331) as 



1 + 



dp- 



fo(P) 



( 



+ 



An dp- 



/„'(/>) 



0. 



The waterbag distribution is now homogeneous in q and symmetric in p 

1 1 



MP) = 



In 2Ap 



(®(p + Ap) - ®{p - Ap)) 



(46) 
(47) 

(48) 
(49) 
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and its derivative is given by 

fo(p) = Y n TKp {6{p + Ap) ~ 6{p ~ Ap)) • (50) 

The second quadratic term in Eq. (1481) vanishes, and one obtains 

0= 1+7T f dp /° (j ' > „ =1 r 7T- (51) 

We finally obtain the complex growth rate 



A = ±^-Api, (52) 

which shows that the waterbag is stable beyond the threshold energy U c = 1/12, since the 
energy of the system is given by U = 

Fig. |3] shows the comparison of this analytical prediction with the numerical results 
obtained for the iV-body simulations of the HMF model: the agreement is excellent. 




Figure 3. Growth rate Re(A) of the instability (full line) as a function of the energy U for the 
HMF model (model ( T27l > with k = 0), as obtained analytically in formula d52l ). The crosses are 
the results of exponential fits of the short-time evolution of the magnetization for the A^-body 
HMF Hamiltonian. 



4. The mean-field (p 4 model 

The second example that we consider is the mean-field ip A model introduced by Desai and 
Zwanzig Il38ll . It is a system where the particles are trapped in an external double- well 
potential, and are in addition coupled via a infinite-range force. It is described by the following 
Hamiltonian 

Wf] = ff dqdpf(q,p) £ + l£-(l-0)£)-l q jj dq'dp' f(q' , p') q' 

Notice that positive (resp. negative) values of the parameter 8 correspond to attractive 
(resp. repulsive) mean-field forces. We have used the same parametrization introduced 
in Ref. |[38l . which can be shown to be minimal by conveniently rescaling the variables 



•(53) 
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and time. The magnetization M is now defined as M[f] = jj dqdp f(q, p) q, so that the 
mean-field potential is given by V[f](q) = ~(9/2)qM[f], whereas the external potential is 
W(q) = q 4 /4-(l -6)q 2 /2. It displays a double well for 9 < 1 and a single well otherwise. The 
solution in the canonical ensemble has been recently derived in Ref. 0411 . emphasizing that 
the system exhibits a second order phase transition. When 9 = 1/2, the critical temperature 
has been found to be T c 0.264, corresponding to a critical energy U* = T c /2 0.132. The 
model has been also solved in the microcanonical ensemble and the entropy as a function 
of energy and magnetization has been derived using large deviations 03 [351 EU, giving 
equivalent results. However, it has been shown that, in the microcanonical ensemble, magnetic 
susceptibility can be negative ll36l [TI. 

For this system, the dispersion relation (|24~1) takes the following form 

-A<f>/tj(J) nq 

f + 0M[f]f (J)——q I e A *' MJ) dq'=0. (54) 
ojkJ) Jo 

The magnetization M[f] can be factored out by multiplying this latter expression by q and by 
integrating it over the phase-space. One gets 

/n e -A<j>lcj(J) nq 
dJfo(J) (J) #-^7p<? 2 J e A *' MJ) dq' = 0. (55) 

Before proceeding to the numerical solution of the above dispersion relation, let us derive 
explicitly the expression that allows us to obtain the stability threshold by setting A = in the 
previous formula. The last integral in Eq. (1551) gives trivially q, while d(p/co can be rewritten 
as dq/p thanks to Eq. (1231) . One finally gets 

1 + 9 J dJf^J) ^dq = 0. (56) 

Let us now restrict to those stationary distributions for which the mean-field vanishes, i.e. 
M[fo] = 0. This case includes those distributions that are symmetric with respect to q = 0. 
For clarity purposes, we shall also restrict to waterbag distributions that have a boundary 
energy U > 0, i.e. fo(J) is constant for all actions < J < J(U) and zero for J > J(U). 
Waterbags with both positive and negative boundary energy U are shown in Fig. [31 
By introducing the following set of variables 

q = xq, (57) 



Ah + (1 - 9) 2 - (1 - 9), (58) 



^4h + a-9) 2 + a-9) 



(4h + (1 - 9) 1 -{1-9) 
the momentum of a particle with positive energy h can be written as 

a 2 

p = + J2(h - W(q)) = V(P 2 -* 2 )0 +x 2 ) . (60) 

V2 

Note that x varies in the range [-p; p], so that the maximum position along a trajectory is pq. 
Now, the action-angles variables (1121131) assume the following form 

J = -2— (£ yj(p 2 -x 2 )(l+x 2 )dx (61) 
2V27T J 



Stability of inhomogeneous states 14 



^3 



V2 



= ^ [(P 2 - 1K(-P 2 ) + (P 2 + D^(-P 2 )] , (62) 

d : =Ai-p\ (63) 

2 ^ J V(P 2 - * 2 )(1 + * 2 ) nq 

(p = ^ T — dx : = to^\-,-p 2 \. (64) 



V(p 2 -x 2 )(l +JC 2 ) \P 
Using the following relation 

; ^ m =y/2q£ , ^ (65) 
V2(/i - J V(P 2 - * 2 )(1 + * 2 ) 

= 4V2^[A-p 2 )-^(-p 2 )], (66) 

one can show, taking also Eqs. (l62l) and (l63~l) into account, that 

; =|dg| = titt 2 — 7T - n — 7T-- (67) 

V2(/i - W(q)) TiP 2 ~ 1) (P 2 - 1)^ 

Considering that g 2 (p 2 - 1) = 2(1 - 0) and q 2 p 2 l(p 2 - 1) = 2/z/(l - 6>), we eventually get the 
following expression for the stability threshold 

Let us now consider the case of the waterbag defined by Eq. (l36l) with J\ = and J 2 = J(U). 
The dispersion relation for this waterbag reads 

1 - 1 3 - 4 ] = 0. (69) 



1 + ^— ^ foW 3/-4--i d/ = 0. (68) 



1 - \ to(U)J(U) 

Solving numerically this latter equation for 6=1/2 gives the threshold energy U c - 0.144, 
which turns out to be pretty close to the value of the statistical transition energy U* found in 
Ref. Il34l . We can also solve numerically the dispersion relation (1551) and obtain the growth 
rate Re(/i). In Fig. |4] this growth rate is compared to a fit of the short-time exponential 
growth of the magnetization obtained by integrating numerically the iV-body Hamiltonian. 
Unfortunately the agreement is only qualitative, although the stability threshold is correctly 
reproduced. 



5. Concluding remarks 

Systems with mean-field interactions are well described by the Vlasov equation in the N —* oo 
limit. An infinity of stationary states exists for such equation and the study of their stability is a 
subject of paramount importance. Many exact results about homogeneous stationary appeared 
in the literature and several stability criteria have been applied. Also, inhomogeneous states 
have been treated, but the study of their stability is more complex Il32l [331 1371 |39l 1401 HTTl . 
Characterizing analytically the stability of stationary solutions of the Vlasov equation will 
have an impact also on the characterization of the slow convergence to equilibrium observed 
in systems with long-range interactions [[Ell2ll3lK|51[6]|, in particular on the study of Quasi- 
Stationary-States (QSS), which are ubiquitous long-lived states in the iV-body dynamics of 
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e u 

Figure 4. Left panel: Representation in the (q, p) plane of the waterbags in action-angle. 
The boundary energies are U = -0.035, 0.05 and 0.2 for lighter and lighter grey levels. The 
dashed line corresponds to the separatrix, which has energy h = 0. Right panel: Growth rate 
Re(A) (full line) computed by numerically solving Eq. OBY The crosses represent the short- 
time exponential rate of growth of the magnetization obtained in numerical simulations of the 
iV-body Hamiltonian with N - 10 6 . 

long-range systems. It has been shown that the lifetime of QSS diverges algebraically with N 
in some simple models and it has been conjectured that this can happen only when the QSS 
corresponds to a stable stationary state of the Vlasov equation HI (see also Ref. Il47ll for an 
interesting mathematical result along this direction). Again, most of the studies on QSS are 
for the homogeneous case. 

In this paper, we have discussed a class of models where, besides the mean-field 
interaction, particles are subjected to an external potential. The effect of the external potential 
is that of creating an inhomogeneity in the spatial distribution. Hence, these models are 
naturally endowed with inhomogeneous stationary states. After rewriting the Vlasov equation 
in action-angle variables, we have shown that some of this inhomogeneous stationary states in 
conjugate coordinates transform into homogeneous stationary states that are homogeneous in 
angle. We have therefore applied the standard tools of linear stability of the Vlasov equation 
to derive a dispersion relation, given in formula (1241) . which is the key result of this paper. We 
have specialized this formula for two models: the HMF model with additional asymmetry ll37l 
and the mean-field <p A model lT34ll35l[36l[38ll . For these two models it is possible to further 
simplify the dispersion relation and to obtain implicit equations that, solved numerically, give 
both the growth rate of the instability and the stability threshold. When the real part of the 
growth rate vanishes, the state is a stable stationary state of the Vlasov equation. We have 
checked these results against the numerical simulation of the Hamiltonian dynamics of the 
corresponding iV-body system. The stability thresholds are in general in good agreement with 
the theoretical predictions, but for the growth rate the agreement is only qualitative for the 
tp 4 model. A case in which the growth rate turns out to be in perfect agreement with the 
simulations is the one of the HMF model |fITl[T2l[T3l 

Those inhomogeneous stationary states that are also stable are good candidates to become 
QSS at finite N. We have therefore pointed out the existence of a new class of inhomogeneous 
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QSS, for which it will be possible in the future to study the law of divergence of the lifetime 
with system size. 
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